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ABSTRACT 

We show that a standard Shakura-Sunyaev accretion disc around a black hole with 
an accretion rate M lower than the critical Eddington limit does show the instability 
in the radiation pressure dominated zone. We obtain this result performing time- 
dependent simulations of accretion disks for a set of values of a and M . In particular 
we always find the occurrence of the collapse of the disc: the instability develops always 
towards a collapsed gas pressure dominated disc and not towards the expansion. This 
result is valid for all initial configurations we tested. We find significant convective 
heat flux that increases the instability development time, but is not strong enough to 
inhibit the disc collapse. A physical explanation of the lack of the expansion phase 
is proposed considering the role of the radial heat advection. Our finding is relevant 
since it excludes the formation of the hot comptonizing corona -often suggested to be 
present- around the central object by the mechanism of the Shakura-Sunyaev insta- 
bility. We also show that, in the ranges of a and M values we simulated, accretion 
disks are crossed by significant amplitude acoustic waves. 
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1 INTRODUCTION 

This work concerns th e possibility of occurrence , stated by 
Shakura and Sunyaev iShakura fc SunvaevlllQT^) . of an in- 
stability in the a-disks when the radiation pressure dom- 
inates, i.e. in the so-called A zone. Shakura and Sunyaev 
demonstrated the existence of thermal instabilities con- 
nected with a difference between Q_ and Q+, i.e. the energy 
emitted per unit area of the disc and the rate of energy- 
generation in the disc. The problem of the existence and 
outcome of the Shakura-Sunyaev instability is important in 
accretion disc physics because it affects the model of the ori- 
gin of the comptonization cloud around some compact ob- 
jects whose spectra contain a significant part in the X-ray 
band. The idea of a hot gas cloud, called 'comptonization 
cloud', around the disc region close to the compact object, 
where photons are pushed towards high frequencies by the 
Compton scattering with electrons, is one of the most com- 
mon ways to explain the spectral behavior in the X band 
(Sunyaev & Titarchuk 1980). In general, the outcome of the 
Shakura-Sunyaev instability is guessed to be the formation 
of a hot cloud around the in ternal disc region, in which 
comptonization could happen iShapiro et alJll976D . 
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Some authors already studied the problem of the a-disc 
time evolution. In 1984 Taam and Lin found that the lo- 
cal Shakura and Sunyaev stability analysis is confirmed 
by global time-dependent simulations of the canonical disc. 
The instability a ppears and gives r ise to luminosity fluctu- 
ations or bursts iTaam fc ^111984 . Taam and Lin's result 
is similar to ours despite of the one-dimensional feature of 
their simulations (whereas we perform two-dimensional ax- 
isymmetric simulations). In 1987 Eggum et al. found that 
an a-disc, at sub-Eddington accretion rates, develops the 
Shakura-Sunyaev instability a nd, as a consequenc e of that, 
collapses in a cold thin sheet jEggum et al.lll987^ ■ In 1998 
Fujita and Okuda simulated a radiation pressure domi- 
nated a-disc in the subcritical acc retion rate regime an d 
found that it is thermally stable llFuiita fc Okudalll99dl . 
They analyzed the disc structure and found that it corre- 
sponds to the configuration of the slim accretion disc model 
(Abramowicz ct al. 1988). In this model, even if the disc 
is radiation pressure and electron-scattering dominated, the 
equilibrium structure is thermally stable, as a consequence 
of an advective radial heat flux that is dominant with respect 
to the vertical radiative heat flux. The cases we present are 
different, since the structure of the disks we simulated is 
close to the Q-model configuration, rather than to the slim 
accretion disc structure. 
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In all the mentioned 2D works the evolution of the system 
has been simulated only for short time-scales, so their au- 
thors cannot exclude they detect only a transient behavior. 
In 2001 Agol et al., by performing local two-dimensional 
hydrodynamic simulations, showed that, for vertically in- 
tegrated dissipation rate proportional to the vertically inte- 
grated total pressure, i.e. in the Shakura- Sunyaev model hy - 
pothesis, the thermal instability develops (lAgol et al.l200l|l . 
The result of this instability is generally the disc collapse. 
Only if a strong enough initial increase of the local radiation 
energy density is given to the disc, there is an expansion, 
that the authors cannot follow because a large amount of 
matter is pushed out of the simulation region. 
In our simulations we confirm the results obtained by Eggum 
and by Agol. Furthermore our studies can follow the system 
evolution on a much larger time-scale and therefore give 
stronger relevance to the collapse result. We can also say 
that the post collapse phase is very long, at least 4.5 10^-^ 
(where Rg is the Schwartzschild gravitational radius of the 
black hole and c is the light speed). We also suggest that 
the convective (in z and r directions) energy transport in 
the disc may have a twofold role: 1) to increase the time- 
scale of the collapse instability; 2) to inhibit the expansion 
instability. 



duj 
vpr — 

dr 



(7) 



V = ctVaH is the kinematic viscosity, a is the viscos- 
ity parameter of the Shakura-Sunyaev model, Va is the local 
sound speed, H is the disc vertical thickness, uj is the local 
angular velocity and the other terms have the usual gas dy- 
namic meaning. 

The gravitational force produced by the black hole is 
given by the well - know n pseudo-newtonian formula by 
IPaczvnski fc Wiital Jl98d) : 



F : 



GM R 

{R-Rg)^R 
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where Rg is the Schwartzschild gravitational radius of 
the black hole, given by: 



Rn 



2GM 



and A4 is the black hole mass. 



(9) 



We adopt the local thermal equilibrium approxima- 
tion for the radiation transfer treatment. However this 
assumption does not affect our conclusions. 



2 THE PHYSICAL MODEL 

The time dependent equations describing the physics of ac- 
cretion disks are well known. They include: 
mass conservation 



Dp 

--— = — p div V 
Dt ' 

radial momentum conservation 
DVr \^ / ,. 

p— — = + 9r+ (dlV (7)r+fr 

vertical momentum equation 

Dv^ _ IdP 
Dt p dz 

energy balance 



-gz + fz 
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Here is the comoving derivative, F is the radiation 
flux, given by: 



F= 



V Erad 



(5) 



3p(fc + (jt) 

k and ctt are the free- free absorption and Thomson scat- 
tering coefficients, Erad is the radiation energy per unit vol 

ume, / is the radiation force, given by: 



7 k + ar 
P F 



(6) 



A is the angular momentum per unit mass, E — e + 
"^p"'' is the total internal energy per unit mass, including 

gas and radiation terms, a is the viscosity stress tensor. 
The component of a that is important in accretion disks is 
the r-(j) one, given by: 



3 THE NUMERICAL METHOD AND THE 
SIMULATIONS PERFORMED 

We set up a new version of the Smoothed Particles Hy- 
drodynamics (SPH) code in cylindrical coordinates, for axis 
symmetric problems. We remind that SPH is a lagrangean 
interpolating method. Recently it has been shown it is equiv- 
alent to finite element s with spars e grid nodes moving along 
the fiuid flow lines For a detailed account 

of the SPH algorithm see .MonagharJ ^1985^. For cylindr i- 
cal coordinates implementati on see [iVIolteni et alJ il99d) : 
IChakrabarti &: Moltenil lll993f) . Our code includes viscosity 
and radiation treatment. Let us note that -in general- a la- 
grangean code is better suited to capture convective mo- 
tions than eulerian codes. With the same spatial accuracy 
(cell size equal to particle size) the SPH particle motion is 
tracked with great accuracy, i.e. the particle size may be 
large but its trajectory can be still determined 'exactly'. To 
integrate the energy equation we adopted the splitting pro- 
cedure. In the LTE condition the radiation energy density 
changes according to the well known diffusion equation given 
by: 



dErad 



= -div _F=V 



SpKtoi 



V Erad 
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where ktot = k + ax- 

In cylindrical coordinates r, z: 



dErad 



dErad 



QErao 



C d 

r dr \ SpKtot dr J ' r dz \ Sputot dz 



(11) 



The SPH- version of the radiation tran sfer term is given 
following the criteria bv lBrookshawl lll994) . The cylindrical 
coordinate version is given by: 



m =iV!!^r^i^U,||i.v,w^. (12) 



\ dt J i Ti ^ rj 



Pj 
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where for clarity we did not put the subscript rad in 
Erad and where: 



0(13) 



This formula can be obtained by the same procedure 
explained by Brookshaw, but taking into account that -in 
cylindrical coordinates- the particles masses are defined as 
mfc — 2-n pk Vk Ark Az^, that explains the further division 
by rj in the term . 

The reference units we use are Rg, for length values, 
and Rg/c for time values. 

We performed several simulations, the ones commented here 
had the following parameter values: 

a) a = 0.01, M = 0.7, domain Ri - R2 = 3 - 300, 
h = 0.3; 

b) Q = 0.1, M = 0.3, domain R1-R2 = 3-200, h = 0.5; 

c) a = 0.01, M = 0.7, domain Ri ~ R2 = 42 - 58, 
h = 0.04; 

where M is in units of Me and Me is the critical 
accretion rate. For all cases the central black hole mass is 
M = 10 Mq The spatial resolution we adopt is h. In the 
case 'a' we have — 29704 particles at time t — 0. 



W e used a variable h procedure iNelson fc Papaloizod 
Il994h . The h values above reported are the initial ones. In 
our procedure, in order to have a not too small particle size 
(and therefore not too great CPU integration times), we 
put a floor for the h values: h is chosen as the maximum be- 
tween the value given by the variable h procedure itself and 
1/10 of the disc vertical half thickness. So we have nearly 10 
particles along the disc half thickness even in the collapsed 
region. This floor was not adopted for the case 'c' since we 
already have a good resolution with a not exceeding number 
of particles (iV = 67518). 

The boundary conditions of the simulations are not 
flxed: as the SPH particles move around, the simulation re- 
gion follows the form assumed by the disc and the values of 
the physical variables at the boundary of the disc are the 
values that characterize the boundary particles at a certain 
time. 

The spatial extension of the initial configuration is decided 
by establishing a radial range of physical interest and a ver- 
tical extension given by the disc thickness of the Shakura- 
Sunyaev model. 

For radiation, the boundary conditions we used are based 
on the assumption of the black-body emission and partic- 
ularly on the Brookshaw approximation ( Brookshaw 1994) . 
At every time step boundary particles are identified by geo- 
metrical criteria (the particle having the maximum absolute 
value of z in a vertical strip of radial width given by h is a 
boundary particle). The boundary particle loses its thermal 
energy according to the formula given by Brookshaw (that is 
an approximation of the diffusion equation at the single par- 
ticle level), that states the particle cooling rate proportional 
to 

QT 

h? 



(14) 







where 
4acT^ 

3pK,tot 



150 
radius 



Figure 1. The r-z profile of the disc of case 'a' is shown at the 
time t = 282 Rg/c . Every SPH particle is represented by a small 
dot. On the x axis the r values in units of Rg are represented. On 
the y axis the z values in units of Rg are represented. 
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10 - 
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(15) 



Figure 2. The r-z profile of the same disc case 'a' at the time 
t = 12000 Rg/c is shown. 



In all our simulations the boundary particles never 
reach an optical thickness lower than 10. 

Let us now comment the results of our study referring 
to the figures of the simulation data. 

Fig. 1 shows the disc structure for the case 'a' at the 
adimensional time t — 282. The disc profile shows the Z 
height of the disc, constant along the A zone. In this figure 
it is evident an initial large convective motion in the disc in 
a region close to the black hole. 

Fig. 2 shows the same disc at the larger time t — 12000. 
The collapsed zone is clearly shown. 

In Fig. 3, the ratio Prad/Pgas at the times t = 282 and 
t = 12000 is shown, exhibiting evidence of the collapse. 

Fig. 4 shows that the disc luminosity is steadily 
decreasing from the initial theoretical value to much 
lower values due to the lower temperature reached by the 
collapsed A zone. 
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150 200 
radius 

Figure 3. The ratio Pradl Pgas at the times t = 282 and 
t = 12000 is shown. On the x axis the r values in units of Rg 
are represented. The configuration at the later time, collapsed, ex- 
hibits a gas pressure dominated zone up to r = 70i?g , whereas the 
other configuration, at a earlier time, shows that, up to 150 Rg, 
the disc is radiation pressure dominated. 




15000 
time 



Figure 4. The time behavior of the disc luminosity is shown. 
On the X axis the time values in units of Rg/c are represented. 
On the y axis the luminosity values in units of erg sec~^ are 
represented. 

Case 'b' has a larger h and a lower number of particles 
and it was possible to follow on the simulation up to the 
large time t = 680000. The r — z distribution of the disc 
particles is very similar to case 'a'. 

We show in Fig. 5 the luminosity of this disc ver- 
sus time. It is apparent that the luminosity has a large 
decrease during the collapse going down from a peak at 
L = 6.5 10^^ erg sec~^ to L = 8.5 10^^ erg sec"\ After 
the collapse the disc luminosity has a slow increase. Nearly 
after the time t = 4.8 10^ (corresponding to physical time 
t = 48 sec, for our parameters) the luminosity reached an 
average value of L = 3 10^^ erg sec~^ and starts to show 
a strong flaring activity with luminosity reaching values up 
to L = 1.6 10^ erg sec~^ . Obviously the possibility of a 
'recharge' of the inner zone is to be expected and has been 





100000 200000 300000 400000 500000 600000 700000 

time 

Figure 5. The time behavior of the disc luminosity is shown for 
the disc of case 'b'. 




Figure 6. The border particles of the disc case 'c' are shown 
at the times i = (vertical crosses), t = 20000 Rg/c (diagonal 
crosses) and t = 22400 Rg/c (asterixes). 



guessed by Eggum and Agol, but -up to now- it was not 
given any estimate of the time-scale involved. It has to be 
noted that the refilling time-scale is definitely shorter than 
the radial viscous drift time-scale tdrift ~ r^/v as derived 
from the canonical disc structure. 

In the 'c' case we simulated a small sector of the disc in 
the radiation pressure dominated zone with a very accurate 
spatial resolution. This case is similar to the Agol's case 1 
(Aeol ct al. 2001), with nearly the same disc parameters and 
numerical resolution (the finest grid used for their case 1 is 
a 256 X 256 one, we have 67518 particles). We also obtained 
a vertical collapse instability of the disc and no expansion. 

Fig. 6 shows the r — z profiles at times t = 20000 and 
t — 22400 of the collapsing disc, compared with the disc 
profile in the initial configuration. The collapsed structure 
is completely gas pressure dominated. In the figure is also 
apparent the fact that there are waves travelling along the 
disc. 
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Figure 7. The radiative, convective and heating time-scales, or- 
dered from the bottom to the top of the figure, are shown versus 
r. On the y axis the time values in units of Rg/c are represented. 

During the simulations, we also calculate the z- 
averaged convective and advective heat fluxes at different 
values of r. With 'convective flux' we mean the heat flux 
along the z-direction due to the vertical motion of the gas, 
whereas with 'advective flux' we mean the heat flux along 
the r-direction due to the radial motion of the gas. The 
main time-scales involved in our problem are defined as 
follows. 

The mean convection time-scale is tconv = , where 

Fconv is the z-averaged convective heat flux and e is the 
total energy density. The mean radiation diffusion time- 
scale is trad = "f^i where Frad is the z-averaged radiative 
heat flux and is the radiation energy density. The mean 
heating time-scale is theat = '^"^1°' , where dErot is the 

rotational energy of the disc ring between r and r + dr, Tr^ 
is the z-averaged viscosity stress and u) is the local angular 
velocity of the disc. 

Fig. 7 shows the mentioned time-scales evaluated at dif- 
ferent radii r = 10 — 100 Rg for case 'b', after time t — 10^. 

From the comparison among the evaluated three time- 
scales it is clear that the convection time-scale is -in general- 
intermediate between the other two time-scales, i.e. the con- 
vection time-scale is greater than the vertical radiation dif- 
fusion time-scale and smaller than the heating time. So we 
can say that the vertical convective heat transfer may have a 
signiflcant role in the vertical redistribution of the generated 
heat. 

The role of the advection (radial transport) appears 
more relevant than the (vertical) convection. The time vari- 
ation of the convective, advective and radiative fluxes is os- 
cillatory and the oscillation frequency is close to the local 
keplerian frequency (Figs. 8 and 9). 

The approximate equality between the flux oscillation 
and keplerian rotation frequencies is explainable if one as- 
sumes that the flux oscillation is due to the propagation 
of an acoustic wave in the disc. Milsom and Taam already 
showed that accr etion disks are crossed by acoustic waves 
l|Milsom fc Taaml fl996. 199't). In the limit of negligible vis- 
cosity and advection and for propagation direction per- 



1e+019 




-1e+019 I ' 1 ' ' 1 1 ' 1 

585000 585500 586000 586500 587000 587500 588000 588500 589000 

time 

Figure 8. The convective and advective fluxes at r = 30 Rg 
are shown versus time. On the y axis the fluxes values in units 
of erg sec~^ cm~^ are represented. The dashed line represents 
the convective flux behavior, whereas the continuous one refers 
to the advective flux. 
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Figure 9. The convective and advective fluxes at r = 80 Rg 
are shown versus time. On the y axis the fluxes values in units 
of erg sec~^ cm~^ are represented. The dashed line represents 
the convective flux behavior, whereas the continuous one refers 
to the advective flux. 



pendicular to the disc rotation axis, it is possible to have 
an analytical formula for the wave frequency. The disper- 
sion relation for these waves is given by = k'^ + K^v^ 
JChandrasekhailll96J) . where ui is the angular frequency of 
the wave, K is the wave number and k is the so-called local 
epicycle frequency, given by k^ = 2^-£:{r^ft), where Q is 
the local angular velocity of the disc, k is in general close to 
the local keplerian angular velocity. Therefore the fact that 
the waves frequencies and the keplerian angular velocities 
are close is a suggestion of the acoustic origin of such waves 
in accretion disks. 

Finally, the waves we obtain have amplitudes that re- 
main roughly constant in time. No other kind of disc struc- 
ture variability is present. 
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We may add that the wave phenomenon appears -at 
different intensity levels- in all cases we examined. We never 
obtained a monotonic regular radial speed as the Shakura- 
Sunyaev model predicts. 



4 DISCUSSION 

Here we discuss two items : the absence of the expansion 
instability and the refilling of the inner zone of the disc. 

The origin of the instability in radiation pressure dom- 
inated zone is clear, but it is not clear why the preferred 
evolution is towards the collapsed state and not towards the 
expanded one. At subcritical accretion rates, i.e. in our ac- 
cretion rate regimes, the expansion instability never occurs. 
We suggest that this result can be due to two effects: the en- 
hanced cooling in the expansion branch and the significant 
role of advection. 

If the disc evolved towards the expansion instability the 
disc density would go down and therefore also its optical 
thickness. The basic LTE approximation, under which the 
disc model is built, breaks down; the disc can then radiate its 
energy content more quickly than it is heated by its viscosity. 
In our view the role of convection and advection may also 
contribute to reduce the expansion instability, but not the 
collapse instability. This mechanism has been considered ca- 
pable to reduce the Shakura-Sunyaev instability. It is known 
that such a transfer should be present in accretion disks 
jRisnovatvi-Knga,n Rlinnikovll 977l:lFuiita Okudall994 
IShakura et alJll97a) . as a result of the entropy gradient as- 
sociated with the radiation field in thermal equilibrium with 
the gas in the accretion disc. This gradient assumes a very 
large value when the temperature distribution is determined 
only by the radiative heat transfer and the viscous energy 
production rate. Following Bisnovatyi-Kogan and Blinnikov 
analysis, the vertical convective motion and the consequent 
vertical heat transfer have, therefore, the role of producing 
an isentropic z structure. If we indicate with Fconv the ver- 
tical convective heat flux, the equation of the vertical heat 
transfer is, in presence of convection; 



^(-Frad + Fconv) + 







(16) 



where Frad is the standard radiative fiux. So we have a 
heat flux whose value has been enhanced from Frad to: 

Fconv 



Ftot — Frad H~ Fconv - — -frad(t H~ ' 



= fFr, 



where 

conv 



(17) 



(18) 



Frad 

It has been demonstrated that this enhancement of 
the vertical heat flux by a factor / produces an in- 
crease of the instability dev elopment time by the same fac- 
tor / JShakuraet a lj (l978ll . From our simulations we ob- 
tain that, in the parameter range we explored, / is al- 
ways less than 2. Such a small value of / does not in- 
crease significantly the instability development time. How- 
ever the advected r adial heat is not taken into accoun t 
in the calculations oflBisnovatvi-Koean fc Blinnikov! il977t) ; 
IShakuraet al]lll978ll . 

Figs. 10 and 11 show the radial and vertical speeds of 
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Figure 10. The radial speed of the disc parcels is shown ver- 
sus r. On the y axis the speed values in adimensional units are 
represented. 
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radius 



Figure 11. The vertical speed of the disc parcels is shown ver- 
sus r. On the y axis the speed values in adimensional units are 
represented. 



the disc parcels versus r. It is clearly apparent from Figs. 8, 
9, 10 and 11 that the radial fiux is greater than the vertical 
one. 

As the radial heat flux assumes a significant value, the 
heat produced at a certain radius is carried away from that 
point of the disc where it had the possibility to produce the 
expansion instability. Furthermore this advection helps the 
cooling of the disc towards the collapsed conflguration of the 
instability. 

We obtain values of the radial heat flux that are large 
with respect to the vertical convective flux. Therefore we 
propose that the absence of the disc expansion may be due 
also to the role of the radial advective heat transfer. 

As regards the collapse instability development times, 
we have compared the theoretical values with the ones 
calculated from the luminosity time behavior: we see that, 
when the disc collapses, its luminosity decreases very 
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5 CONCLUSIONS 




Figure 12. The ratio Prad/ Pgas in the initial configuration (as- 
terixes) and at the times t = 10^ (diagonal crosses) and t = 6 10^ 
(vertical crosses) is shown. The initial configuration, in the shown 
radial range, is radiation pressure dominated. At the intermediate 
time (t = 10^) the disc, collapsed, is gas pressure dominated. The 
configuration at the later time {t = Q 10^), refilled, shows again 
a radiation pressure dominated zone up to r = 60 Rg . 



quickly and the time on which the luminosity reduction 
occurs can be assumed as the instability development time. 
As an example, we can see the comparison between the 
two times in the case 'b'. The theoretical instability de- 
velopment time, tth, can be calculated from the instability 
development rate f2t?i, given by Q,th = qcj , with 

A{l3r) = 8 + 51/3r — 3/3^ and ui equal to the a ngular velocity 
of th e disc at the radius considered (.Shakura fc SunvaevI 
Il976l) . With this expression for Q.th we can calculate the 
instability development time tth as ttu ~ jf-^- The value of 
tth at the innermost radii is of the order of 10~^ sec. 
From the luminosity time behavior we deduce similar 
values. 



Fig. 12 shows the ratio of the radiation pressure to the 
gas pressure versus the radial distance for the case 'b' when 
the disc was in the initial configuration, in the collapsed state 
(time t = I 10^) and when it was in the refilled state (time 
t = 6 10*). It is clear that, in the A zone, the refilled state has 
a Prad/ Pgas valuc larger than 1, but smaller than the very 
initial value, at time t — 0. The presence of this zone with 
Prad > Pgas also iu the refilled configuration suggests that 
the Shakura-Sunyaev instability may be again operating to 
produce the recurrent fiaring activity. 

The refilling time-scale computed from the simulation 
data is about 5 10^ , that is shorter than the theoretical drift 
time-scale due to the viscosity. We argue that the refilling is 
driven not only by the viscosity: also the wave phenomenon 
may play significant role. Indeed the plain theoretical model 
predicts a time-scale tdrift = jv = 1.8 lO'^. This fact 
should be taken into account developing theoretical models 
of recurrent flaring disc activity due to refilling. 



The results of our simulations of the time evolution of aPtot 
disks, with a large portion in radiation pressure dominat- 
ing conditions, show definitely that the Shakura-Sunyaev 
collapse instability develops. After a long time the disc re- 
covers partially its luminosity and furthermore it exhibits a 
flaring like activity. The recovery time-scale is shorter than 
the viscous drift time-scale. We attribute to the advective - 
convective heat transfer a signiflcant role to determine both 
the recovery and the not occurence of the expansion insta- 
bility. 

We suggest that the reason for which we see no expan- 
sion instability is the presence of a radial convective motion, 
that isn't considered in the Shakura and Sunyaev analysis, 
but is naturally simulated by our code. In such conditions, 
the radial heat transfer due to the advection of matter (and 
thermal energy) carries the excess of energy produced by vis- 
cosity (that would cause the thermal expansion) away from 
the disc element at that r. These results imply some problem 
for the model of the formation of the comptonization cloud 
supposed to exist in many disc configurations to explain 
the observed radiation spectra. From the canonical disc it 
seems impossible to produce that comptonization cloud via 
the plain Shakura-Sunyaev instability. 

Another relevant aspect of our simulations is the pres- 
ence of acoustic waves; they can cause a kind of periodical 
variability of the radiation spectrum and intensity that could 
be related, from the ob servational p oint of view, with the 
phenomenon of QPOs JPsaltislliool . Oscillation frequen- 
cies of the radiative flux we calculated are compatible with 
some QPOs. It is true that the total disc luminosity takes 
into account the contributions from the whole disc and so 
the wave phenomenon may not appear. However one should 
take also into account the possibility that the reflection of 
the radiation coming from inner disc zones by outer zone 
waves could produce an enhanced oscillation in the total lu- 
minosity. Such a study requires an accurate treatment of the 
interaction between the emitted radiation and the disc itself 
and of the radiation transfer in the outermost layers of the 
disc, that is beyond this contribution. 

The flaring activity could be also responsible of QPO 
phenomena. We are planning to confirm and study its char- 
acteristics by numerical simulations with increased spatial 
resolution. 
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Figure captions 

Fig. 1. The r-z profile of the disc of case 'a' is shown at 
the time t = 282 Rg/c . Every SPH particle is represented 
by a small dot. On the x axis the r values in units of Rg 
arc represented. On the y axis the z values in units of Rg 
are represented. 

Fig. 2. The r-z profile of the same disc case 'a' at the time 
t = 12000 Rg/c is shown. 

Fig. 3. The ratio Prad/Pgas at the times t = 282 and 
t — 12000 is shown. On the x axis the r values in units of 
Rg are represented. The configuration at the later time, 
collapsed, exhibits a gas pressure dominated zone up to 

r = 70 Rg, whereas the other configuration, at a earlier 
time, shows tiiat, up to 150 Rg, the disc is radiation 
pressure dominated. 

Fig. 4. The time behavior of the disc luminosity is shown. 
On the X axis the time values in units of Rg/c are rep- 
resented. On the y axis the luminosity values in units of 

erg sec~^ arc represented. 

Fig. 5. The time behavior of the disc luminosity is shown 
for the disc of case 'b'. 

Fig. 6. The border particles of the disc case 'c' are shown at 

the times t = (vertical crosses), t = 20000 Rg/c (diagonal 

crosses) and t — 22400 Rg/c (asterixes). 

Fig. 7. The radiative, convective and heating time-scales, 

ordered from the bottom to the top of the figure, are shown 

versus r. On the y axis the time values in units of Rg/c are 

represented. 

Fig. 8. The convective and advective fluxes at r = 30 Rg 
are shown versus time. On the y axis the fluxes values 
in units of erg sec~^ cm~^ are represented. The dashed 
line represents the convective flux behavior, whereas the 
contirmous one refers to the advective flux. 
Fig. 9. The convective and advective fluxes at r = 80 Rg 
are shown versus time. On the y axis the fluxes values 
in units of erg sec^^ cm^^ are represented. The dashed 
line represents the convective flux beiiavior, whereas the 
continuous one refers to the advective flux. 
Fig. 10. The radial speed of the disc parcels is shown versus 
r. On the y axis the speed values in adimensional units are 
represented. 

Fig. 11. The vertical speed of the disc parcels is shown 
versus r. On the y axis the speed values in adimensional 
units are represented. 

Fig. 12. The ratio Prad/Pgas in the initial configuration 

(asterixes) and at the times t — 10^ (diagonal crosses) and 
t = 6 10^ (vertical crosses) is shown. The initial config- 
uration, in the shown radial range, is radiation pressure 
dominated. At the intermediate time {t = 10^) the disc, 
collapsed, is gas pressure dominated. The configuration at 
the later time {t = 6 10^), refilled, shows again a radiation 
pressure dominated zone up to r = 60 Rg. 
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